pacman::p_load(tidyverse, estimatr)

df <- read_csv('datasets/ypccc.csv')
m1 <- lm_robust(generaterenewable~sea_level_rise + density + pct_college +
                  pop_18_35 + pop_35_64 + pop_65_plus + 
            pct_white + inc_20_50 + inc_50_100 + inc_over_100 +
            pct_dem, df)
m2 <- lm_robust(prioritycleanenergy~sea_level_rise + density + pct_college +
                  pop_18_35 + pop_35_64 + pop_65_plus + 
            pct_white + inc_20_50 + inc_50_100 + inc_over_100 +
            pct_dem, df)
m3 <- lm_robust(localofficials~sea_level_rise + density + pct_college +
                  pop_18_35 + pop_35_64 + pop_65_plus + 
            pct_white + inc_20_50 + inc_50_100 + inc_over_100 +
            pct_dem, df)
m4 <- lm_robust(president~sea_level_rise + density + pct_college +
                  pop_18_35 + pop_35_64 + pop_65_plus + 
            pct_white + inc_20_50 + inc_50_100 + inc_over_100 +
            pct_dem, df)
m5 <- lm_robust(prioritycleanenergy~sea_level_rise + density + pct_college +
                  pop_18_35 + pop_35_64 + pop_65_plus + 
            pct_white + inc_20_50 + inc_50_100 + inc_over_100 +
            pct_dem, df)
m6 <- lm_robust(rebates~sea_level_rise + density + pct_college +
                  pop_18_35 + pop_35_64 + pop_65_plus + 
            pct_white + inc_20_50 + inc_50_100 + inc_over_100 +
            pct_dem, df)
m7 <- lm_robust(reducetax~sea_level_rise + density + pct_college +
                  pop_18_35 + pop_35_64 + pop_65_plus + 
            pct_white + inc_20_50 + inc_50_100 + inc_over_100 +
            pct_dem, df)
m8 <- lm_robust(regulate~sea_level_rise + density + pct_college +
                  pop_18_35 + pop_35_64 + pop_65_plus + 
            pct_white + inc_20_50 + inc_50_100 + inc_over_100 +
            pct_dem, df)
m9 <- lm_robust(fundrenewables~sea_level_rise + density + pct_college +
            pop_18_35 + pop_35_64 + pop_65_plus + 
            pct_white + inc_20_50 + inc_50_100 + inc_over_100 +
            pct_dem, df)

extractCoef <- function(x){
  return(x %>%
    tidy %>%
    filter(term == 'sea_level_rise'))
}

dat <- bind_rows(
  extractCoef(m1),
  extractCoef(m2),
  extractCoef(m3),
  extractCoef(m4),
  extractCoef(m5),
  extractCoef(m6),
  extractCoef(m7),
  extractCoef(m8),
  extractCoef(m9)
) 

metadf = 
  data.frame(
    effects = dat$estimate,
    ses = dat$std.error
  )
out <- meta::metagen(TE = metadf$effects, 
        seTE = metadf$ses, 
        sm = "SMD", 
        random = TRUE, 
        common = FALSE, 
        method.random.ci = "HK")

dat <- bind_rows(
  dat,
  data.frame(
  estimate = out$TE.random,
  conf.low = out$lower.random,
  conf.high = out$upper.random,
  outcome = 'Pooled Effect'
))

dat %>%
  mutate(outcome = factor(outcome, levels=c(
    'Pooled Effect',
    'fundrenewables',
    'regulate','reducetax','rebates',
    'generaterenewable','prioritycleanenergy',
    'localofficials','president'
  )),
  z = factor(c(0,0,0,0,0,0,
        0,0,0,1))) %>%
  ggplot(aes(x=outcome, y=estimate, ymin=conf.low, ymax=conf.high, shape=z)) +
  geom_pointrange(size=0.5,fill='white') +
  geom_hline(yintercept=0, linetype=2, color='red') +
  coord_flip() + theme_bw() +
  labs(x='',y='Coefficient')  +
  scale_shape_manual(values = c(21,24)) +
  theme(legend.position='none')
ggsave(width=4, height=2,filename='figures/ypccc_regression_results.pdf')
